be a good place for these?

Yeah, the more the marrier.

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/test/interface/mass_matrix_tests.jl is probably a good test to slam some around.

Turns out that lots of these are still not inferrable--even with this PR, going from Rosenbrock DAE AD tests, the move from

sol = solve(prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8)

to

using ADTypes: AutoForwardDiff
sol = @inferred solve(prob_mm, Rodas5P(autodiff=AutoForwardDiff(chunksize=3)), reltol = 1e-8, abstol = 1e-8)

already runs into an issue, since the type of solve is still inferred to be Any. I think there is something else going on here (or I don't know how to use @inferred properly, which is possible), but I can't tell what makes this different from the MWE I've been testing against.

That said, even though this PR apparently doesn't do enough for the existing test suite, it is already helping with the MWE and therefore probably my code base.

It's fine to sprinkle some @test_brokens around to get things in, and document them.

OK, so the quirk in the test I was looking at turns out to be that for

M = Diagonal([1.0, 1.0, 0.0])
roberf = ODEFunction(rober, mass_matrix = M)

get_differential_vars(roberf, u0) is fully inferrable, since M was constructed as diagonal. But in the test, we had

M = [1.0 0 0
     0 1.0 0
     0 0 0]
roberf = ODEFunction(rober, mass_matrix = M)

With a full matrix, inference's best guess at get_differential_vars(roberf, u0) is Union{BitVector, DifferentialVarsUndefined}.

We currently get DifferentialVarsUndefined if all of the following are satisfied:

• Matrix is not a UniformScaling
• Matrix has at least one zero element
• Matrix is either an AbstractSciMLOperator or not diagonal

So, for example, a matrix like [1 1 0; 0 1 0; 0 0 1] with one off-diagonal nonzero currently gives DifferentialVarsUndefined.
What should be the exact criteria for having the differential vars undefined? My instinct is that even a fully-zero matrix should still just give us a falses(size(mass_matrix, 1)), so maybe this branch should not even exist?

What should be the exact criteria for having the differential vars undefined? My instinct is that even a fully-zero matrix should still just give us a falses(size(mass_matrix, 1)), so maybe this branch should not even exist?

I believe the proper definition is that the number of algebraic variables is determined by the size of the null space of M. If M is not diagonal the algebraic variables are the e_i vectors which live in the null space? I think you can always orthogonalize a basis of the null space down to just being defined by linear combinations the algebraic variables.

But I couldn't figure out how to calculate this easily so the diagonal cases were handled and any case where someone couldn't figure out the differential variables but needed it turned into an error. In lots of cases it's not necessary to calculate so it ended up being an okay solution, with a few special cases like dense. But yes... I don't know, I think in general you probably need to QR factorize the M and then take the Q part and check its null space?

julia> M = [1 1 0; 1 1 0; 0 0 0]
3×3 Matrix{Int64}:
 1  1  0
 1  1  0
 0  0  0

julia> qr(M).R
3×3 Matrix{Float64}:
 -1.41421  -1.41421  0.0
  0.0       0.0      0.0
  0.0       0.0      0.0

That definition isn't quite right either.

But I couldn't figure out how to calculate this easily so the diagonal cases were handled and any case where someone couldn't figure out the differential variables but needed it turned into an error. In lots of cases it's not necessary to calculate so it ended up being an okay solution, with a few special cases like dense.

OK, that makes sense--I had the sense my solution might be too easy, otherwise it would have already been there. I'll revert it to keep the previous behavior, I think

For at least some of the tests, I can just change to explicitly using a Diagonal matrix, and do likewise in the examples to point people this way, plus maybe add a note on this type-stability question to the mass-matrix docs. That would be easy enough, and that way we avoid doing any possibly-expensive linear algebra (i.e., for cases with larger mass matrices) in what should typically be a simple preparation step.